% GAUTIER LE BIHAN - 2020
% Replication files for "Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector Menu-Cost Model?" Review of Economics and Statistics
%
% Figure 5 + Figure M CalvoPlus
clear;
tic
cd ..\..\Simulations_VC\Counterfactual_sim\MS_produits_CalvoPlus
load ..\..\Simulations_VC\MS_produits_CPlus\stat_d
test=stat_d;

load ..\..\Estimation_param\MS_produits_CPlus\actual_moments_k
weight=actual_moments_k(:,12);
sum(weight)
secteur=actual_moments_k(:,1);
test_act=[actual_moments_k(:,2:6)];


% Suppose lambda= cte
load ..\..\Simulations_VC\Counterfactual_sim\MS_produits_CalvoPlus\stat_simu_out.mat

for jj=1:3
    
    param=stat_simu_out(jj).stat_simu_out(:,5:9);
 end 
secteur=actual_moments_k(:,1);

f=param(:,1);
f_p=param(:,2);
med=param(:,3);
interq=param(:,4);
kur=param(:,5);

moments_sim=[f f_p med interq kur];
moments_actual=[test(:,5:10)];



 figure(1);
%P0
subplot(2,4,1);

%param=stat_simu_out(1).stat_simu_out(:,5:9);
moments_sim=[test(:,5:9)];
moments_actual=[test_act(:,1:4)];
x=moments_sim(:,1);
y=moments_actual(:,1);
sz=weight*150;
sz=20
scatter(y, x,sz, 'o', 'k');
mdl = fitlm(x,y)
mat_R2(1,1)= mdl.Rsquared.Adjusted

x1=[0  0.2]
y1=x1;
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1]);
ylabel('Freq. sim.') % x-axis label
xlabel('Freq. data') % y-axis label
title('Baseline')

subplot(2,4,2);

param=stat_simu_out(1).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,1);
y=moments_actual(:,1);
sz=weight*150;
sz=20
scatter(y, x,sz, 'x', 'k');
mdl = fitlm(x,y)
mat_R2(1,2)= mdl.Rsquared.Adjusted

x1=[0  0.2]
y1=x1;
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1]);
ylabel('Freq. sim.') % x-axis label
xlabel('Freq. data') % y-axis label
title('\lambda')

subplot(2,4,3);
param=stat_simu_out(2).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,1);
y=moments_actual(:,1);
sz=weight*150;

sz=20
scatter(y, x,sz, '+', 'k');
mdl = fitlm(x,y)
mat_R2(1,3)= mdl.Rsquared.Adjusted

x1=[0  0.2]
y1=x1;
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1]);
ylabel('Freq. sim.') % x-axis label
xlabel('Freq. data') % y-axis label
title('\mu')

subplot(2,4,4);
param=stat_simu_out(3).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,1);
y=moments_actual(:,1);
sz=weight*150;
sz=20
scatter(y, x,sz, '*', 'k');
mdl = fitlm(x,y)
mat_R2(1,4)= mdl.Rsquared.Adjusted


x1=[0  0.2]
y1=x1;
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off


axis([x1 x1]);
ylabel('Freq. sim.') % x-axis label
xlabel('Freq. data') % y-axis label
title('\sigma')



subplot(2,4,5);
moments_sim=[test(:,5:9)];
moments_actual=[test_act(:,1:4)];
x=moments_sim(:,3);
y=moments_actual(:,3);
sz=weight*150;

sz=20
scatter(y, x,sz, 'o', 'k');
mdl = fitlm(x,y)
mat_R2(2,1)= mdl.Rsquared.Adjusted

x1=[0 0.15]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('Median sim.') % x-axis label
xlabel('Median data') % y-axis label
title('Baseline')


subplot(2,4,6);


param=stat_simu_out(1).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,3);
y=moments_actual(:,3);
sz=weight*150;

sz=20
scatter(y, x,sz, 'x', 'k');
mdl = fitlm(x,y)
mat_R2(2,2)= mdl.Rsquared.Adjusted

x1=[0 0.15]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('Median sim.') % x-axis label
xlabel('Median data') % y-axis label

title('\lambda')

subplot(2,4,7);

param=stat_simu_out(2).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,3);
y=moments_actual(:,3);
sz=weight*150;

sz=20
scatter(y, x,sz, '+', 'k');
mdl = fitlm(x,y)
mat_R2(2,3)= mdl.Rsquared.Adjusted

x1=[0 0.15]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('Median sim.') % x-axis label
xlabel('Median data') % y-axis label
title('\mu')
subplot(2,4,8);

param=stat_simu_out(3).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,3);
y=moments_actual(:,3);
sz=weight*150;
sz=20
scatter(y, x,sz, '*', 'k');
mdl = fitlm(x,y)
mat_R2(2,4)= mdl.Rsquared.Adjusted

% hold on
% xd=y;
% scatter(y,xd,'.','k')
% hold off
x1=[0 0.15]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off

axis([x1 x1])
ylabel('Median sim.') % x-axis label
xlabel('Median data') % y-axis label
title('\sigma')



 print('..\..\..\figures\figure7_ref1.pdf','-dpdf', '-fillpage')
 print('..\..\..\figures\figure7_ref1','-depsc')

 
 
 
 
 
 
 
 
  
 figure(2);
%P0
subplot(3,4,1);

%param=stat_simu_out(1).stat_simu_out(:,5:9);
moments_sim=[test(:,5:9)];
moments_actual=[test_act(:,1:4)];
x=moments_sim(:,2);
y=moments_actual(:,2);
sz=weight*150;
sz=20
scatter(y, x,sz, 'o', 'g');
mdl = fitlm(x,y)
mat_R2(3,1)= mdl.Rsquared.Adjusted

x1=[0.5 1]
y1=x1;
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1]);
ylabel('Fracup. sim.') % x-axis label
xlabel('Fracup. data') % y-axis label
title('Baseline')

subplot(3,4,2);

param=stat_simu_out(1).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,2);
y=moments_actual(:,2);
sz=weight*150;
sz=20
scatter(y, x,sz, 'o', 'k');
mdl = fitlm(x,y)
mat_R2(3,2)= mdl.Rsquared.Adjusted

x1=[0.5  1]
y1=x1;
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1]);
ylabel('Fracup. sim.') % x-axis label
xlabel('Fracup. data') % y-axis label
title('\lambda')

subplot(3,4,3);
param=stat_simu_out(2).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,2);
y=moments_actual(:,2);
sz=weight*150;

sz=20
scatter(y, x,sz, 'o', 'b');
mdl = fitlm(x,y)
mat_R2(3,3)= mdl.Rsquared.Adjusted

x1=[0.5  1]
y1=x1;
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1]);
ylabel('Fracup. sim.') % x-axis label
xlabel('Fracup. data') % y-axis label
title('\mu')

subplot(3,4,4);
param=stat_simu_out(3).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,2);
y=moments_actual(:,2);
sz=weight*150;
sz=20
scatter(y, x,sz, 'o', 'r');

mdl = fitlm(x,y)
mat_R2(3,4)= mdl.Rsquared.Adjusted

x1=[0.5  1]
y1=x1;
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off


axis([x1 x1]);
ylabel('Fracup sim.') % x-axis label
xlabel('Fracup data') % y-axis label
title('\sigma')



subplot(3,4,5);
moments_sim=[test(:,5:9)];
moments_actual=[test_act(:,1:4)];
x=moments_sim(:,4);
y=moments_actual(:,4);
sz=weight*150;

sz=20
scatter(y, x,sz, 'o', 'g');
mdl = fitlm(x,y)
mat_R2(4,1)= mdl.Rsquared.Adjusted

x1=[0 0.15]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('IQR sim.') % x-axis label
xlabel('IQR data') % y-axis label
title('Baseline')


subplot(3,4,6);


param=stat_simu_out(1).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,4);
y=moments_actual(:,4);
sz=weight*150;

sz=20
scatter(y, x,sz, 'o', 'k');
mdl = fitlm(x,y)
mat_R2(4,2)= mdl.Rsquared.Adjusted

x1=[0 0.15]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('IQR sim.') % x-axis label
xlabel('IQR data') % y-axis label

title('\lambda')

subplot(3,4,7);

param=stat_simu_out(2).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,4);
y=moments_actual(:,4);
sz=weight*150;

sz=20
scatter(y, x,sz, 'o', 'b');
mdl = fitlm(x,y)
mat_R2(4,3)= mdl.Rsquared.Adjusted

x1=[0 0.15]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('IQR sim.') % x-axis label
xlabel('IQR data') % y-axis label
title('\mu')
subplot(3,4,8);

param=stat_simu_out(3).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
x=moments_sim(:,4);
y=moments_actual(:,4);
sz=weight*150;
sz=20
scatter(y, x,sz, 'o', 'r');
mdl = fitlm(x,y)
mat_R2(4,4)= mdl.Rsquared.Adjusted

% hold on
% xd=y;
% scatter(y,xd,'.','k')
% hold off
x1=[0 0.15]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off

axis([x1 x1])
ylabel('IQR sim.') % x-axis label
xlabel('IQR data') % y-axis label
title('\sigma')


subplot(3,4,9);
moments_sim=[test(:,5:9)];
moments_actual=[test_act(:,1:5)];
x=moments_sim(:,5);
y=moments_actual(:,5);
sz=weight*150;

sz=20
scatter(y, x,sz, 'o', 'g');
mdl = fitlm(x,y)
mat_R2(5,1)= mdl.Rsquared.Adjusted

x1=[0 7]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('Kurto sim.') % x-axis label
xlabel('Kurto data') % y-axis label
title('Baseline')


subplot(3,4,10);


param=stat_simu_out(1).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
moments_actual=[test_act(:,1:5)];

x=moments_sim(:,5);
y=moments_actual(:,5);

sz=weight*150;

sz=20
scatter(y, x,sz, 'o', 'k');
mdl = fitlm(x,y)
mat_R2(5,2)= mdl.Rsquared.Adjusted

x1=[0 7]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('Kurto sim.') % x-axis label
xlabel('Kurto data') % y-axis label

title('\lambda')

subplot(3,4,11);

param=stat_simu_out(2).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
moments_actual=[test_act(:,1:5)];

x=moments_sim(:,5);
y=moments_actual(:,5);

sz=weight*150;

sz=20
scatter(y, x,sz, 'o', 'b');
mdl = fitlm(x,y)
mat_R2(5,3)= mdl.Rsquared.Adjusted

x1=[0 7]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off
axis([x1 x1])
ylabel('Kurto sim.') % x-axis label
xlabel('Kurto data') % y-axis label
title('\mu')
subplot(3,4,12);

param=stat_simu_out(3).stat_simu_out(:,5:9);
moments_sim=[param];
moments_actual=[test(:,5:10)];
moments_actual=[test_act(:,1:5)];

x=moments_sim(:,5);
y=moments_actual(:,5);
sz=weight*150;
sz=20
scatter(y, x,sz, 'o', 'r');
mdl = fitlm(x,y)
mat_R2(5,4)= mdl.Rsquared.Adjusted

% hold on
% xd=y;
% scatter(y,xd,'.','k')
% hold off
x1=[0 7]
y1=x1
hold on 
plot(x1,y1, 'Color',[0.5 0.5 0.5])
hold off

axis([x1 x1])
ylabel('Kurto sim.') % x-axis label
xlabel('Kurto data') % y-axis label
title('\sigma')
 print('..\..\..\figures\figure7append_ref1.pdf','-dpdf', '-fillpage')
 print('..\..\..\figures\figure7append_ref1','-depsc')

 